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UNSTEADY SPHERICAL FLOW BEHIND A KNOWN SHOCK LINE 


ABSTRACT 


‘ The hydrodynamical equations of unsteady spherical flow are 
converted into characteristic form and solved numerically by a 
difference method. The "initial-value" curve is the shock line 
obtained by the least-square fit to some compiled shock-front data 

on spherical Pentolite, of such form as to approach Kirkwood-Brinkley's 
theoretical asymptotic shock-front decay curve. Results are tabulated 
on positive sound paths, mass particle paths, and lines of constant 
distance. 
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Superscript: 
i 
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Subscript: 
° 
1,2,3,66. 


Roman 


hits. 


SYMBOLS 


represents iteration index 


denotes dimensional quantity 


denotes state of the undisturbed air 
identify points in the t,m space 


radius of charge 
local sound velocity = %(p%, s* J 
o* 5* 


c* 
ct 


specific internal energy 
e* - et 
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p* 


specific enthalpy = e* + 5 
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* 
% 
a tou” 


4& function proportional to the mass between a particle 


“path and the path of the boundary between alr and c.- 


_ plosion gas (see eq. 11-7) 


re 


Rta 


total pressure 
ye. 
p "Py 
v5 
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radial distance 


el 
at 


gas constant 


: s* = specific entropy . 
4 s* - s* 
a 8 = 2 
' RF 
‘ s. refers to the shock line in the t,m space 
m - 
t* = time 
% 
t= — t* 
4 Te = absolute temperature 
a [* 
T = a 
% 
ut a mass velocity of air 
ue 
us ow 
oO * 
UH = shock velocity 
ye ‘ 
Uc ~— 
% 
Woe . | 
har’pudt ; constant m 
] “shock 
Greek; - 
Om refers to the forvard~facing sound path in the t,m@ space 
Oy ar - refers to the forvard-facing sound path io the t,r space 
om au refers to the forvard-facing sound path in the p,u space 
| Bem refers to the backward-facing sound path in the t,m space 
Ey / : 
| 4 Per refers to the backward-facing sound path in the t,r space 
A Bau refers to the backvard-facing sound path in the p,u space. 
. 
“f . 
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"t,m refers to the particle path in the t,m space 


tyr refers to the particle path in the t,r space 
4 %,u refers to the particle path in the p,u space 
: p* = density 
q % 
Po 
% Oe 
Ww = 1.4 oe oF 
coMe) 


I. INTRODUCTION 


Extensive experimental data, obtained both by measurement of shock- 

: front arrival times and by piezo-electric gage measurement of hydrostatic 

: pressure change across the front, exist for the propagation of the blast 
wave from bare sphericgl charges exploded in air. These data are sufficient 
Yor the algebraic determinat‘on of all the remaining air flow parameters 

| such as density, entropy, and particle velocity across the shock front. 

| Parameters behind the front, however, cannot be measured as accurately nor 

| be calculated as easily. The pressure-time curve recorded by the piezo- 

Hi electric gage has scatter in the data in hoth the pressure p and the time t* 

directions. But if these observed pressure-time records are assumed 

correct and known at each distance ry i.e., if the function 


p = Be, r) 
is assumed known, then the particle velocity u’= u(¢,2) and the entropy 


a= s( tr) can be calculated by integratica of a cet of ordinary differential 
5 equations (Ref. 1). 


Here, we consider the numerical integration for the air flow parameters 
pu, and é, given the shock line only, in the cone of determinacy of 
this line. 


The system of partial differential equations deseribing this flow is 
of the hyperbolic type. At cach point in the ft, r spaee there exist 


three characteristics, two of which are path lines of sonie disturbances 
travelling outward und inward relative to the fluid particles, and the third 
of which is the path line of particles. This nature of the differential 
équations limits the domain of determinacy of the shock line; in particular, 
this domain cannot be extended into the gone of explosion gases without 
additional information about the zone. 


The method of calculation consists of replacing the original system 
of fiow equations by characteristics, which are written in Pinite difference 
form and solved numerically by step-by-step construction of the characteristic 
netvork. Computations are performed in the BRL electronic high-speed 
computer GORDVAC. . 


To tee 


THECEDINO BAO BLUMK.NOT Mua 


ia 


il 


eee, 


a eee Ae 


. ain. 


The shock-line chosen as the “initial curve" is an analytic expression 
fitted to some compiled experimental data on uncased spherical Pentolite 
charges fired in air at standard conditions away from reflecting obstacles. 
The smoothness of this curve necessarily destroys some information about 
discontinuities in the flow field. 


The results are tabulated as functions of time along particle paths, 
positive sound paths, and constant radius lines. Some representative 
results are given graphically. These results may be not only of practical 
interest in the analysis of blast measuring apparatus and in the study 
of damage to structures, but may also be of theoretical interest in 
approximating blast-wave calculations. 


II. EQUATIONS OF FLOW 


The differential equations describing the continuous spherical flow 
of a non-conductive and non-viscous compressible fluid are, in Eulerian 
coordinates (Ref. 2), 


Conservation of mass: 


(II-1) se + wu a + p* ae + ae = 0; 


Conservation of momentum: 


II-2) du* ou* 1 opt ; 
( sx“ * u* + Sa a QO 5 


ore. Ge 
Adiabaticity: 
(II-3) ds* ds* 
ot WO Sw eo: 


The independent coordinates are time t* and distance r*; the dependent are 
pressure p*, mass velocity u*, specific entropy s*, density p*, and sound 
velocity c*. ‘These equations are supplemented by the equations of state 
(11-4) ce = ct (pe, ut) , 

pt = pt (pt, at) , 
which are tabulated for alr dn Ref. 3. By the total derivative of the 
equations above 

dp* eo a dp* e spc ipes ct) dg® 

Z oe 8 


and by (II-3), we may put (II-1) into the form 
9 
3 ; : pe hye 
(1-5) oe u* SES + cope a + co ao. 
In non-dimensional form equations (11-5), (II-2) and (11-3)tecome, respectively, 


Se * wu at * we Su + Seu % Oo. 
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(11-6). ¢ du ou ce 3 
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The Lagrangian form of these equations is obtained by replacing 
the independent variable r by a variable related to the mass bounded 
within the radial coordinate of each particle. We define a variable 
m such that 


a 
or e 

om __ aru 

ot C ° 


The compatibility condition oe = an for continuous second derivatives 
reduces to tne conservation of mass equation, and is automatically satisfied. 
m, so defined, is proportional to the fluid mass between two spherical 
shells of radii ry and y moving with the fluid. In particular, we let ry 
be the radius of the boundary between the explosion gas and air. ‘The 


transformation (II-7) uow converts the flow equations (II-6) into 


or 2.2 au @ucu 
Toes Soe eS es 


au 2 dp 

(11-8) Sire 
hae 
xs” 


in the t, m space. 
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He TII. CHARACTERISTIC EQUATIONS 

Be The system of partial differential equations (II-6) or (II-8) is of 


the hyperbolic type; i.e., three real characteristics exist in the t,r 
or t,m space along which discontinuities in derivatives can propagate 
(Ref. 2). We denote these characteristics by a, 6B and y, with subscripts 
to specify the space. The characteristic equations then are: 


a ar. 
P iG a at =u + c . 
aa (III-1) _ dm 9 2 
.: Sm ae OF 
= (o4 : L op, gy acu, 
EF Pw w at at a 
a ar 
aa Bort ae 7 Uo. 
a dm _ 2 
(III-2) Boat ap eae 
=e ee . 2 ap au 2enu , 
. ane put a a ees i, 
ae : ar 
a = "tyr at uy 
ae - - : Gm _ 0 ' 
ae * “(XII-3) ‘t,m' at ’ 
f ee ds 
i : oem 5S 0 » 
a3 . "pai at 
7 ote _ Physically, the q@ and 6 characteristics correspond to the forward 
Be ey and: backward facing sound paths respectively, and y. corresponds to the 
a ae particle path. ‘Solving this system of characteristic equations is equi- 


Pea Oso 2. valent ‘to tolving the original set of flow equations, (II-8) or (II-6). 
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IV. SHOCK-FRONT CONDITIONS AND INITIAL DATA 


Across the discontinuous shock front travelling into stationary air 
with velocity U*, the parameters of flow are related by the following 
Rankine-Hugoniot equations of conservation of mass, momentum and energy 
(Ref. 2). 


Conservation of mass: 
p¥(UK - uk) = pk UK; 


Conservation of momentum: 
px(UF - ut)*+ pA = ox ure 4 Px ; 


Conservation of energy: 


ye 
1. 2 ede f) 
g (ue - ut)" + et + Be 5 Us * eo * oF ° 


These equations are derivable from (II-6) as weak solutions. We supplement 
the above equations by the equations of state 


p¥ = p#(pt, at) , 
e* = e&(pt,s*) , 
c¥ = c#(pt, at) , 


where e* is the specific internal energy. Shear and Day have tahulated 
these equations for air in Ref. 43. In the undisturbed state 


a us. af, the air is assumed to obey the ideal gas law ps s Rtpore 
with specific heat ratio 1.4. 


In dimensionless form the equations abecve, together with the equations 
cf the shock front path, may be written 


Dn _ om _ Om Dr | 2 
De ot > thee U; 
Dr 

p - UV» 

= (U + u) =zl4u, 


(3v-1) 
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e=e(p,s) , 


where D/Dt is differentiation along the shock line. The nine unknowns 

in the eight independent equations in (IV-1) above are m, r, p, u, 8, U, 

®, ¢, @. Thus, an experimental observation of any of these variables as 

a function of t or r determines the remaining variables as functions of 
t or Yr. 


Extensive experimentation has been conducted on the propagation of 
shock waves from spherical Pentolite, because of the reproducibility of 
ite explosion characteristics. Two methods of measurement have been 
commonly employed: in the first, chock arrival times between several 
points of observation are observed as functions of distance in the form 
t= t(r), ty either photographic observation of shock front path in the 
t,r plane or by piezo-electric gage observation of sudden pressure changes; 
in the second method, magnitudes of pressure jump across the shock are | 
measured by means of plezo-electric gages with calibrated voltage output. 
Goodman (Ref. 4) has compiled these data and constructed the empirical fit - 


r-l_ 2 ri 1}* + hy 
(1v-2) p = M8.a6 (5-2 ro + oh7.0 7 


which approaches Kirkwood-Brinkley's asymptotic solution (Ref. 5) 


f = (rin ny“ ‘for r -» 6 (A oni Bare constants). The fitting is based 
on an abundance of points fur 1< r <€ 100, and ona relatively small 
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number for 100<r<200. Kirkwood-Brinkley's asymptotic curve with para- 
meters determined at r = 150 differs at most about 10 percent from this 
expression out to our maximum distance r = 8200. We somewhat arbitrarily 
ccnsider r : 150 to be the demarkation point beyond which the shock Line 
de equivaleny to Kirkwocd-Brinkley's asymptotic solution. Hence, the 
domain of determinacy of the experimental data lies between the positive 
characteristics @= 0 and @ = 91 passing respectively thru the initial 

a ir = 15° points on the shock line (Fig. 1). 


From the empirical equation above and the shock conditions (IV-1) 
and also from the Hugoniot table of Ref. 3 are derived the complete shock-~- 


_ line values. We represent these values functionally by 


m= mt’, 
r= r(t), 
(Iv-3) 8: p= pt), 
us u(t) ,) 
po a(t), 


which we assume to be continuous and to possess continuous first cerivatives. 
This set of functions will be regarded as the "inicdal-value" line for the 


numerical integration of the hydrodynamical aquavions of flow. 
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V. METHODS OF NUMERICAL INTEGRATION 


e) 
Since, for shock velocity U>0O, the slopes 1.4 r°U of Sem? 
a 


re 4 2 2 +, 
4 or” of Oem -wr” of By om and O of "ym? are respectively related by 


q 1 or > La rt, 0> - are 

(Ref. 2), no characteristics are tangent to the shock line; the domain 
of determinacy of the shock line in the t,m space is therefore the area 
bounded by the shock line Sy om » the particle path y = 0 through the 
initial point of ae » and the forward-facing sound path m= 156 
through the terminus of Se om (Fig. 1). 

With the specified “initial data" we solve the flow equations by 
converting the equivalent characteristic equations into finite difference 
ia form. Though methods using differences higher than the first can be 
q : employed, machine limitations restricted the calculation here to first 
; difference methods only in which the lattices are kept small. The 
coefficients of the differences are averaged, and improved after each 
3 : cycle of the iterative process. The final results at each point are 
; : then accurate to the third order in lattice size, as we see later. 


The air torence scheme used is specifically as follows. In Figure 2, 
let 1B and 3 be newents oF adjacent aa lines given by the characteristic 
equations (IIT-1), th und 33 be segments of ned jacent Vig lines given 
by the characteristic equations (i11-4), und 45 a ego of a fi, 
line given by the characteristic equations (ITI-2). 
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Figure 2. Construction of the Characteristic Lattice. 


Of these five points let 1, 2 and 5 be completely known, let 4 be the 
point to be calculated, and let 5 be known to the (i - 1)'th iterative 
cycle. We denote the iterative cycle by superscript (1) and location 

of the points by numerical subscripts. The i'th values at point 4 are 
then calculated by means of the difference forms of the sets of character- 
datic equations (1ZI-1), (111-2) and (I1I-3). The coordinates (1), al!) 
are calculated from the difference equations 


al) on (wr*), G2) * (ar), 


tm 2 er 


an - 5 
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(1-1) 


The values of tS are now compared. If a - t,2 0 , 


and te 


(1) (4) , ; 5% and 45 
then te ,m," are calculated as the intersection of 23 and 45: 


ae nt) : n(*) - (ax? )(t2) i (ax®) (4-4) 
tym i I Shee Ca eg ee a 
ten > ty 
(v-2) 
(1) 


Tent FR 
If al - t,< 0, then 42) ‘ al?) are calculated as the intersection 


om “7 
of 12 and 45: 


ous al) . n(*) (wr?) (4-2) f (wr® (4-2) 
a Saye 
aati" 
(v-3) On! nf?) = n(e{!) ; 


where nf?) 5 n(t{?)) ie either an interpolation formula or the differ- 
ence equation cig 13. That is, if point 5 is on 23 in the (4 - 1)'th 
iteration cycle, the "t, ‘a curve through point 2 is used to find the 

i'th value of point 5; but if the (1 - 1)'th value of point 5 is on 2, 
the a. curve through point 2 is used. All other required quantities 
at point 5 are obtained either by the interpolation formulas 


yg v(+f*)) 
(1) xe) 
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(v-4) ug) a u(ef)) 


aft) fehl!) 
(4) o(f!)) 
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“~*~ 7~ 
along 23 or 12, or by the difference equations along these arcs, 


depending on the location of point 5. (Machine limitation required 
a) 

that we represent 12 and 23 by linear interpolation formulas in our 

computation. ) 


An alternative procedure for calculating point 5 is to use either 
the set of 6 and y equations only or the set of B and @ equations 
only, rather than to alternate between the two sets depending on the 
previous iterative value of t53 but in this case the accuracy of the 
constructed Beem curve through point 4 ae decreases with increasing 
distance of point 5 from the lattice Tosh. Another methOr for cal- 
culating point 5 is to use the @ equations and the curve 3 repre- 
sented by interpolation formulas between points 1 and 3; but in this 
method machine :ound-off difficulties will be encountered in regions 
where the geometry of the lattice Tash is such that the curve if is 
nearly tangent to either of the characteristics Th or 3h. 


With point 5 known, p and u at point 4 can be calculated from 


(3) (4-2) ’ (3 ) pi?) 


a a ns oe oo 
pu 
: 1 ts - ts 


(ea fe + (&S )s 


ae a~ 


| Gis GY eft) - a 
Bp. 4 My 
p,u ee a 


pe ee eee ee ee Nr ok Kieth r= pany he OLS p i ne) 5s 


w and c at point 4 are obtained from the equations of state 


oe aol), al! 
(v-7) 
of) = ofp,{), aft) 
Thus, the i 'th values of all the variables oft) , a?) : (4), 
pf), ul), of?) and oft ) are determined. The process is repeated 


until the (i + 1)'th Re the i'th values differ insignificantly. 


In the above procedure, when i = 1 the O'th iterative values at 
point 4 required for averaging of coefficients of the difference 
equations can be obtained either by extrapolating along the computed 
a. ore through point 3 down to a m or by averaging the coefficients 
at points 1 and 3. For example, we can set, as is done in this report, 


(ar®), + (ur*) 


2 9) 
(ax®)() = (®)() 2» — 2, 
For the se that we require in the comparison of a) with to 


vefore point 5 16 calculated, we can similarly use the average 
2f0) 3 Pree See - : ° 


A special case of the above procedure for calculating point 4& 
occurs in the neighborhood of the shock line. Here, the procedure 
is the same as at interior points, except that points 2 and 4 are 
considered ecineident (Fig. 3), and interpolations (V-4) along TB are 
earried out along the shock Line 3. i rather than along ‘em OF Be thru 2. 
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Figure 3. Characteristic Lattice Adjoining the Shock Line 


The points integrated by the preceding scheme may be indexed by 
assigning successive values to the members of the @ and 7 characteristics, 
such that on the shock line these indices are equal. The characteristic 


‘network in the @, 7 space would then be as given in Fig. 4. The known 


points on the shock line S are (1,1),°(2,2), (4,5), ete., and the 
unknown points below S are calculated in the sequence (1,0), (2,1), 


(2,0), (3,2), (5,1), (3,0), ete. 


(156 , 156) 


1,0) (1,0) (2,0) (30) (1560) o—» 


Pigure 4. Lattice ina , 7 Space 


Once the lattice poinus along the shock line are chosen, the 
lattice network in the domain of integration becomes determined, unless 
the structure of the network io arbitrarily changed somewhere in the 
domain. Therefore, to minimize round-off errors, the shock-line lattice 
size in the t,m space is wade variable. In the decaying shock con- 


_oidered here, the lattice interval At along the shock line io made an 


increasing function of t. 


The integration process described uses Te ny and a. a equations 
é ’ d 


to determine the characteristic network points in the t,m space; t.e., 


the lattices are rectangular in the a, 7 space in regions away fron 

the shock line. Thus, integration results, including values on fae 
bomndaries of the domain of determinacy, becom: tabulated on @ and y 
characteristics, and values on f characteristics must be obtained by 
interpolation of the results. If, novever, results alorg other combinations 
of characteristics are preferred, the process can be inodified. We can, 

for example, construct the vharactertstic network in thea, 6 space. 

In Fig. 2, if pointy 5 and 3 are known, differeace forms of the Ge in 

and Be in characteristic equations iv (I1I-b) aad (I1I-5) respectively 

can be used to determine the quantities ts My. In the special case 
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that values along y are not of interest, this method has the advantage that 

only 3) need be interpolated as a function of mM,» whereas all interpolations 
(V-4) are required in the previous method; however, it has the disadvantage 

of requiring a special procedure at the boundary y = O separating air 


Gi ies a A aes act one 
a 


from the explosion gas. 


A still another method is to choose the integration points arbitrarily 
along constant t lines, say t = t, + At (Figure 5), and to employ the 
intersection points 1, 2, and 3 of the Oe om? By om? em characteristics 
respectively with the known line t = t to calculate quantities at the 
integration point 4. This procedure is permissible because t = constant 
lines are nowhere tangent to any characteristic in any continuous and 
finite domains, However, several disadvantages arise. Firstly, the 
lack of automatic sdjustment of network size according to gradients, as 
dn the methods that construct the characteristic network, necessitates 


the determination of Am along each constant tt line according to the 


-_ results obtained for previous times, in order to control errors. Secondly, 
: ioterpolations must be performed at all three points 1, 2 and 3, whereas 
: previous methods required this operation at one point only. Thirdly, a 


special procedure is required to determine the 7 = 0 boundary passing 
through the initial point of the shock line. Finally, values along 
Cheracteriatics that may be desired for theoretical studies must all be 


subsequently interpolated for. 


Figure 5. Integration Along Constant t Lines. 
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The existence theorem shows that integration methods based on 
differencing the original system of flow equations (II-8) rather than 
differencing the characteristic system require that the unknown point 
be within the domain of determinacy of the known points used in the 
difference scheme (Ref. 5). Since this condition must be examined at 
every integration point, complete use of the characteristic method is 
simple: and probably more accurate where the integration zone is 
fairly smooth, such as occurs in this report. 


FN SN NET SN NT I LITE © Aggy ae oh Ty EL ATE DEST eam Se BS ee ae 8 ae re on 


VI. INTEGRATION ERRORS 


The sets of characteristic equations (III-13, (IIi-2) and (III-3) 


, are in the form 


(vI-1) Sgt) alt) 6) 3609) (xt),(2) (6) 
. Dae f = (¢ peeegp ) gt (g peoesp ); 


jel 


where gi) represents one of the set of six variables (t, m, r, p, u, 8), 
alk represents one of the three characteristic variables (a, 6, 7), 
and superscript 4 refers to one of the spaces (t,r), (t,m) or (p,u). 


=f we assume 


pl dke ) gik4 ) and god) to be expandable in convergent series in the arc 
vWvW-_o 
13 along the @ . characteristic, we have 
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: (32) 
ato = pt) + Taleo) + aT ( xR 2 (2a) )? tree, 


(jk4 seel dks ) 


BA). of Bed), dee pleat) + B- “eee (al)? 5.0. , 


(kL) 
Gt). fe) 2, Ge agit! Yoo") «5, Ce, (aa)? 


Bo 


ght). itt) b ; ar, (aa) + 5 ar sy (aa)? 
() gd MO) ody 2 ae) (x) 2 
x" = by tr ye (40 ) + Che (ar) 


2.(3) 


(3) 
(8) 2 gO nce (aa"*)) + d Ci. (2a'*))? 


from which we obtain by addition and subtraction 
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kt), g(t) okt) , 1 Ci, (ta!*))? + .., 


(3) - gd) G) 33 ( 
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4 
ae olka) , aT ; Oe ) (aq)? 


. jh 45) 
Solving these equations for 2° i) : git and ay), 


and substituting into (VI-1), we obtain 


: Gh) Cie) ka) (ke) 
(vi-2) oe | ae a oe ro : ot (aa"*) 
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This result indicates that for lattices sufficiently small the method of 
first-order differencing with average coefficients applied to the set of 
characteristic equations introduces a truncation error of third order in 
lattice size. While we can further reduce this error by integrating 
with several lattices and extrapolating the results to zero lattice 
(Ref. 6), the experimental accuracy of the shock line does not make 

the process worthwhile. 


The relative round-off error iz kept small by increasing the time 
step At along the shock line in accordance with the diminishing gradients 
of the variables. 


The numerical study of error propagation and growth in scattered 
limited regions of the domain of integration by introduction of small 
changes in the variables, indicates that relative errors remain approx- 
imately of constant order. 


VII. DISCUSSION OF RESULTS 


The comparison between Brode's calculation on TNT (Ref. 7) based 
on the method of artificial viscosity and our results on Pentolite 
calculated by the method of characteristics, must be semi-qualitative, 
since the explosives and the nature of the initial and boundary values 


differ, and since Brode's scaling factors for time and distance depend 


iS 
i 


on uncertain explosion energies ¢*. Brode's dimensionless time and 
distance are defined by 


% 
— re : re 
“Oo Ps 


If \ = 0.0156 is the charge surface r = 1, as Fig. 16 in Brode's report 


i seems to indicate, we can relate 1,A to t, r by 


2. 4 
iGome a. hoes A 


Brode's curves transformed to t,r space are compared with our assumed 


r. shock line and computed data in Fig. 6. (Because of some inconsistencies 
a between Brode's curves in Ref. 7, we have used only Fig. 16-21 in this 

ij reference, supplemented by some of his unpublished data on isobaric lines 
i that he has kindly furnished us.) 


As Fig. 6 indicates, the calculated results influenced by the 
asymptotic segment of the shock line are questionable. While the error 
may be due to fitting the asymptotic curve to inaccurate shock line 
data near the limit of observation, wore probably the asymptotic formula 
itself requires modification in the direction of fastex shock decay, bo 
that pressure at any fixed distance as a function of time decays faster 
in accord with observation, and the expanding contact surface reverses 
its direction of motion, as it eventually must. 


Typical graphs are drawn for the various variables along lines of 
constant @, y, and r. 


The tables give the shock line used as initial values, and the results 
along lines of constunt a, 7, and r in the domain influenced by tne exporie 
montal shock line onty. 
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